Códigos

Ejemplos de cópulas

library(copula)
# cópula independiente
persp(indepCopula(), pCopula, theta = -30, phi = 20)

# Comonotonicidad
n.grid <- 26
u <- seq(0,1,length.out=n.grid)
grid <- expand.grid("u[1]"= u,"u[2]"= u)
M <- function(u) apply(u,1,min) #Cota superior M
x.M <- cbind(grid,"M(u[1],u[2])" = M(grid)) #Evalua M en el grid
wireframe2(x.M)

# Contramonotonicidad
n.grid <- 26
u <- seq(0,1,length.out=n.grid)
grid <- expand.grid("u[1]"= u,"u[2]"= u)
W <- function(u) pmax(0,rowSums(u)-1) #cota inferior W
x.W <- cbind(grid,"W(u[1],u[2])" = W(grid)) #Evalua W en el grid
wireframe2(x.W)

Cotas para cópulas

norm.cop <- normalCopula(0.4)
par(mar = c(4, 4, .1, .1), mfrow = c(1,2))
persp(norm.cop, pCopula, cex.axis = 0.5)
persp(norm.cop, dCopula, cex.axis = 0.5)

Familias arquimedianas

Cópulas de Frank

Cópula de Frank

copfr <- function(x,th){-log((exp(-th*x)-1)/(exp(-th)-1))}
par(pty="s")
curve(copfr(x,th = 5),      
      from = 0.1, to = 1,
      ylab = expression(phi[Fr](x)), col = 1)
curve(copfr(x,th = 1),
      from = 0.1, to = 1, add = T, col = 2)
curve(copfr(x,th = 0.01),
      from = 0.1, to = 1, add = T, col = 3)
curve(copfr(x, th = -1),
      from = 0.1, to = 1, add = T, col = 4)
curve(copfr(x, th = -5), 
      from = 0.1, to = 1, add = T, col = 5)
legend("bottomleft", 
       lty = rep(1, 5),   
       paste("th=", c(5, 1, 0.01,-1,-5)), col = 1:5)

Simulación de la cópula de Frank

library(copula)
set.seed(1)
par(mfrow = c(3,3))
theta <- c(-90, -50, -10, -1, 0 , 5, 20, 50, 90)
for(i in 1:9){
  U <- rCopula(n = 200,
    copula = archmCopula(family = "frank", param = theta[i]))
    plot(U,
         xlab = expression(u[1]), 
         ylab = expression(u[2]),
           main = eval(substitute(expression(paste(theta, " = ", j)), list(j = as.character(theta[i])))))
}

Cópula de Clayton

Función generadora

Phi_C <- function(u, theta = 0.0001)(u^(-theta)-1)/theta
unitario <- seq(0, 1, by = 0.1)
    
curve(Phi_C, from = 0.1, to = 1, col = "red")
for(theta in seq(0, 10, length = 100)){
            curve(Phi_C(x, theta = theta), 
                  from = 0.001, to = 1, add = T, 
                  col = ifelse(theta > 1, "blue", "orange"),
                  ylab = expression(Phi[C](x)))
}

Simulación de cópula de Clayton

library(copula)
set.seed(1)
par(mfrow = c(3,3))
theta <- c(-.95, -0.5, -0.1, 0.1, 1 , 5, 20, 50, 90)
for(i in 1:9){
    U <- rCopula(n = 200,
                copula = archmCopula(family = "clayton", 
                                     param = theta[i]))
    plot(U, 
         xlab = expression(u[1]), 
         ylab = expression(u[2]),
         main = eval(substitute(expression(paste(theta, " = ", j)), list(j = as.character(theta[i])))))
}

Cópula de Gumbel

Función generadora

Phi_G <- function(u,theta=1)(-log(u))^theta
            
curve(Phi_G, from = 0.1, to = 1, col = "red")
for(theta in seq(1, 20, length = 100)){
        curve(Phi_G(x, theta = theta), 
              from = 0.001, 
              to = 1, 
              add = T, 
              col = ifelse(theta > 10,"blue","orange"))
}

Simulación de la cópula de Gumbel

library(copula)
set.seed(1)
par(mfrow = c(2,3))
theta <- c(1, 1.5, 2, 4, 8, 50)
for(i in 1:6){
    U <- rCopula(n = 200,
                       copula = archmCopula(family = "gumbel",
                                            param = theta[i]))
    plot(U, 
         xlab = expression(u[1]), 
         ylab = expression(u[2]),
         main = eval(substitute(expression(paste(theta, " = ", j)), list(j = as.character(theta[i])))))
}

Ejemplos

Primer ejemplo independencia

set.seed(1)
n <- 1000
sigma <- 0.5 # asumida
#Matriz de covarianzas
Sigma <- sigma^2 * diag(2)
Sigma

     [,1] [,2]
[1,] 0.25 0.00
[2,] 0.00 0.25

#Genera las muestras de dos normales
z <- matrix(rnorm(2*n),nrow=2)
#Transforma para escalar las variables
Y <- sigma*diag(2) %*% z 
X <- exp(Y) #X tiene distribución lognormal
dim(X)

[1]    2 1000

X[1:2,1:3]

         [,1]      [,2]      [,3]
[1,] 0.731084 0.6584845 1.1791029
[2,] 1.096169 2.2202957 0.6634948

par(pty = "s") #gráfico cuadrado
plot(X[1,], X[2,], xlim=c(0,5), ylim=c(0,5), pch=16, cex = 0.5, 
     main = "Lognormales con Independencia")

Segundo ejemplo

rho = 0.7 #correlación
Sigma <- sigma^2*matrix(c(1,rho,rho,1),nrow=2)
Sigma

      [,1]  [,2]
[1,] 0.250 0.175
[2,] 0.175 0.250

# Obten la matriz raíz cuadrada B de la matriz 
# definida positiva Sigma
e <- eigen(Sigma)
v <- e$vectors
B <- v %*% diag(sqrt(e$values)) %*% t(v)
z <- matrix(rnorm(2*n),nrow=2)
Y <-  B %*% z #Transforma para escalar las variables
X <- exp(Y)
par(pty = "s") #haz el gráfico cuadrado
        plot(X[1,], X[2,], xlim=c(0,5), ylim=c(0,5), pch=16, cex=0.5,
        main="Lognormales con Correlación=0.7")

Generación de variables aleatorias usando cópula \(t\)

library(mvtnorm) #para generar t multivariada
set.seed(3); TT <- t(rmvt(n = 1000, sigma = Sigma, df = 1)) #genera una dist. t(1)
U  <- pt(TT, df = 1) #distribución t(1) para uniformes
X <- rbind(qgamma(U[1,],2,1), qt(U[2,],2))
par(mfcol = c(1,3)); plot(X[1,], X[2,], pch = 16, cex = 0.5, main = "Copula")
hist(X[1,], main = "Gamma(3,1)", breaks = 30, prob = T)
hist(X[2,], main = "t(2)", breaks = 100, prob = T)

cor(X[1,], X[2,])

[1] 0.4727318

Simulación de cópulas bivariadas

theta <- 5
c1 <- cbind(runif(1000),runif(1000)) # genera u y v
A <- 1+theta*(1-2*c1[,1])
B <- sqrt(A^2-4*(A-1)*c1[,2])
FGM <- cbind(c1[,1],2*c1[,2]/(A+B))
par(pty="s"); plot(FGM[,1],FGM[,2],pch=16,cex=0.3,main=paste("Muestra de FMG",theta)) 

Paquete copula

Clase cópula

library(copula)
copula.normal4 <- ellipCopula(family = "normal", 
                              dim = 4, dispstr = "un",
                              param = c(0.4,0.5,0.2,0,0.3,0.8))
copula.normal4 #objeto de clase normalCopula

Normal copula, dim. d = 4 
Dimension:  4 
Parameters:
  rho.1   = 0.4
  rho.2   = 0.5
  rho.3   = 0.2
  rho.4   = 0.0
  rho.5   = 0.3
  rho.6   = 0.8
dispstr:  un 

u <- rCopula(200, copula.normal4) #genera observaciones de la cópula construida
cor(u)

          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.3999356 0.5583336 0.2587525
[2,] 0.3999356 1.0000000 0.0091090 0.3456752
[3,] 0.5583336 0.0091090 1.0000000 0.7516383
[4,] 0.2587525 0.3456752 0.7516383 1.0000000

pairs(u, pch = 16, cex = 0.5)

Ejemplo 2

micopula.t3 <- ellipCopula(family = "t", 
                           dim = 3, 
                           dispstr = "toep",
                           param = c(0.8,0.5), df = 8)
micopula.t3 #objeto de clase tCopula

t-copula, dim. d = 3 
Dimension:  3 
Parameters:
  rho.1   = 0.8
  rho.2   = 0.5
  df      = 8.0
dispstr:  toep 

rCopula(5, micopula.t3) #genera cinco observaciones de la cópula construida

          [,1]       [,2]      [,3]
[1,] 0.7561208 0.83477990 0.3862287
[2,] 0.1214745 0.15426384 0.3950312
[3,] 0.1372568 0.06109462 0.1017845
[4,] 0.5601066 0.55269459 0.5608740
[5,] 0.4282554 0.43444576 0.6929481

Ejemplo 3

clayton2 <- archmCopula(family = "clayton", 
                        dim = 2, param = 2)
clayton2 #el programa llama alpha al parámetro

Clayton copula, dim. d = 2 
Dimension:  2 
Parameters:
  alpha   = 2

# Generemos una muestra de ésta cópula:
y <- rCopula(1000,clayton2)

Curvas de nivel

par(mfrow = c(1,2))
contour(clayton2,dCopula) #gráfica de curvas de nivel
plot(y, cex = 0.3)

Ejemplo 4

copula.Frank5 <- archmCopula(family = "frank", dim = 3, param = 5)
micopula <- mvdc(copula = copula.Frank5, 
                 margins = c("norm","pois","gamma"),
                 paramMargins = list(list(mean=10,sd=2),
                                     list(lambda=3),
                                     list(shape=2,scale=4)))
u <- rMvdc(300,micopula) #muestra aleatoria
par(mar = c(1,1,1,1))
pairs(u,pch=16,cex=0.5)

Funciones de distribución y densidad para cópulas

(u <- rMvdc(5, micopula)) # Puntos del dominio

          [,1] [,2]      [,3]
[1,] 12.217118    5 23.825533
[2,] 10.959193    6  5.464355
[3,]  9.983242    3  4.927729
[4,]  9.062085    2  2.944202
[5,]  7.397641    1  2.956918

dMvdc(u, micopula) # puntos de la densidad

[1] 0.0004065068 0.0003107008 0.0055394034 0.0062570787 0.0049991691

pMvdc(u,micopula) # puntos de la distribución

[1] 0.81171345 0.36484532 0.26675940 0.10375430 0.03015304

(u <- rCopula(5,copula.Frank5)) # puntos del dominio

          [,1]      [,2]      [,3]
[1,] 0.7872706 0.8028639 0.8070379
[2,] 0.8792277 0.8008430 0.9959447
[3,] 0.3188677 0.6734840 0.6816883
[4,] 0.4410392 0.2512305 0.4597243
[5,] 0.1723916 0.5241862 0.2036805

dCopula(u,copula.Frank5) # puntos de la densidad

[1] 4.3678780 5.6718894 0.6929495 1.7070607 1.3589326

pCopula(u,copula.Frank5) # puntos de la distribución

[1] 0.63703314 0.74680497 0.28060158 0.17364825 0.08518855

Realización de la distribución conjunta

library(scatterplot3d)
par(mfrow = c(1,2),
    mar = c(1,2,1,1),
    oma = c(0,0,1,1),
    mgp = c(2,1,0))
u <- rMvdc(200, micopula)
scatterplot3d(u, cex.symbols = 0.5, pch = 16)
v <- rCopula(200, copula.Frank5)
scatterplot3d(v, cex.symbols = 0.5, pch = 16)

Contornos de funciones de densidad

miMvd1 <- mvdc(copula = archmCopula(family = "clayton",
                                    param = 2), 
               margins = c("norm", "gamma"),
               paramMargins = list(list(mean = 0, sd = 1),
                                   list(shape = 1, scale = 2)))
miMvd2 <- mvdc(copula = archmCopula(family = "frank",
                                    param = 5.763), 
               margins = c("norm", "gamma"),
               paramMargins = list(list(mean = 0, sd = 1),
                                   list(shape = 1, scale = 2)))
miMvd3 <- mvdc(copula = archmCopula(family = "gumbel",
                                    param = 2), 
               margins = c("norm", "gamma"),
               paramMargins = list(list(mean = 0, sd = 1),
                                   list(shape = 1, scale = 2)))

# Parámetros para tau de Kendall de las tres dist. = 5
par(mfrow = c(1,3), 
    mar = c(2,2,1,1), 
    oma = c(1,1,0,0), 
    mgp = c(2,1,0))
contour(miMvd1, dMvdc, xlim = c(-3,3), ylim = c(0,4))
contour(miMvd2, dMvdc, xlim = c(-3,3), ylim = c(0,4))
contour(miMvd3, dMvdc, xlim = c(-3,3), ylim = c(0,4))

contour(miMvd1, dMvdc,xlim=c(-3,3), ylim=c(0,4))
contour(miMvd2, dMvdc,xlim=c(-3,3), ylim=c(0,4))
contour(miMvd3, dMvdc,xlim=c(-3,3), ylim=c(0,4))

# función persp es similar
persp(miMvd1, dMvdc, xlim = c(-3,3), ylim = c(0,4))
persp(miMvd2, dMvdc, xlim = c(-3,3), ylim = c(0,4))
persp(miMvd3, dMvdc, xlim = c(-3,3), ylim = c(0,4))

Correlaciones

# 1
muestra <- cbind(c(2,3,1,5,4,9,6,4,2,10),
                 c(3,4,5,2,8,6,8,3,1,10))
cor(muestra, method = "kendall")

          [,1]      [,2]
[1,] 1.0000000 0.3953488
[2,] 0.3953488 1.0000000

# La correlación de Pearson usual
cor(muestra)

          [,1]      [,2]
[1,] 1.0000000 0.6440133
[2,] 0.6440133 1.0000000

# 2
X <- c(1,2,3,4,5,6,7)
Y <- c(1,3,6,2,7,4,5)

# 3
muestra <- cbind(c(2,3,1,5,4,9,6,4,2,10),
                 c(3,4,5,2,8,6,8,3,1,10))
cor(muestra, method = "spearman")

          [,1]      [,2]
[1,] 1.0000000 0.5613497
[2,] 0.5613497 1.0000000

#Para efectos comparativos
cor(muestra, method = "kendall")

          [,1]      [,2]
[1,] 1.0000000 0.3953488
[2,] 0.3953488 1.0000000

# La correlación de Pearson usual
cor(muestra)

          [,1]      [,2]
[1,] 1.0000000 0.6440133
[2,] 0.6440133 1.0000000

Ejemplo 3

X <- rexp(100, rate = 3)
Y <- X^3
plot(X,Y)

cor(X,Y)

[1] 0.8494689

cor(X,Y, method = "spearman")

[1] 1

cor(X,Y, method = "kendall")

[1] 1

Aplicaciones

library(copula)
# Define el objeto cópula a generar, en este caso es una normal con correlaciones
# dadas
# El argumento dispstr se refiere a la estructura a la matriz de covarianza que caracteriza a
# la cópula. "un" es para indicar que no tiene estructura.
# ver detalle en https://www.jstatsoft.org/index.php/jss/article/view/v021i04/v21i04.pdf
copula_normal_3 <- normalCopula(c(sin(0.7*pi/2), sin(0.4*pi/2), sin(0.3*pi/2)), dim = 3, dispstr = "un")
set.seed(100) #fija una semilla
U <- rCopula(1000, copula_normal_3) #Genera la muestra aleatoria
pairs(U, pch = 16, cex = 0.5)

round(cor(U, method = "kendall"), 2)

     [,1] [,2] [,3]
[1,] 1.00 0.69 0.39
[2,] 0.69 1.00 0.28
[3,] 0.39 0.28 1.00

W <- cbind(qnorm(U[,1], mean = 4, sd = 5),
qt(U[,2], 4),
qbinom(p = U[,3], size = 25, prob = 0.4))
head(W)

         [,1]        [,2] [,3]
[1,] 2.180664 -0.16821757    9
[2,] 8.355016  0.65959400   11
[3,] 2.275934  0.17518414    8
[4,] 2.902724 -0.09615615   10
[5,] 5.235607  0.58978484   10
[6,] 3.620093 -0.27198890   11

#Grafica los histogramas y agrega densidades con las distribuciones deseadas para ver
#la aproximación
par(mfrow = c(1,3))
hist(W[,1], prob = T, breaks = 50); 
points(sort(W[,1]), dnorm(sort(W[,1]), 4, 5), 
       type = "l", col = "red")
hist(W[,2], prob = T, breaks = 50)
points(sort(W[,2]), dt(sort(W[,2]),4), 
       type = "l", col = "red")
hist(W[,3], prob = T);
points(sort(W[,3]),dbinom(sort(W[,3]), size = 25, prob = 0.4), 
       type = "l", col = "red", lwd = 3)

cor(W, method = "kendall")

          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6884324 0.4120571
[2,] 0.6884324 1.0000000 0.3006960
[3,] 0.4120571 0.3006960 1.0000000

Estimación de una cópula

library(Ecdat) # fuente de datos
library(copula)
library(fGarch)
library(MASS) # usa las funciones fitdistr y kde2d
library(fCopulae) # funciones adicionales de copula (pempiricalCopula y ellipticalCopulaFit)
data(CRSPday, package = "Ecdat")
head(as.data.frame(CRSPday)) # muestra la estructura de los datos

  year month day        ge       ibm     mobil      crsp
1 1989     1   3 -0.016760  0.000000 -0.002747 -0.007619
2 1989     1   4  0.017045  0.005128  0.005510  0.013016
3 1989     1   5 -0.002793 -0.002041  0.005479  0.002815
4 1989     1   6  0.000000 -0.006135  0.002725  0.003064
5 1989     1   9  0.000000  0.004115  0.005435  0.001633
6 1989     1  10 -0.005602 -0.007172  0.008108 -0.001991

ibm <- CRSPday[,5]; crsp <- CRSPday[,7]
(n <- length(ibm)) #número de observaciones

[1] 2528

par(pty = "s")
plot(ibm, crsp, cex = 0.4, pch = 16)
abline(h = 0, v = 0, col="red")

ibm <- CRSPday[,5]
crsp <- CRSPday[,7]
est.ibm <- as.numeric(fitdistr(ibm,"t")$estimate) #parámetros t: media, escala, gl
est.crsp <- as.numeric(fitdistr(crsp,"t")$estimate)
#Convierte los parámetros de escala a desviaciones estándar en el caso de la t
est.ibm[2] <- est.ibm[2]*sqrt(est.ibm[3]/(est.ibm[3]-2))
est.crsp[2] <- est.crsp[2]*sqrt(est.crsp[3]/(est.crsp[3]-2))
#Grados de libertad para cada caso
est.ibm[3]

[1] 4.276156

est.crsp[3]

[1] 3.473982

(tau <- cor(ibm,crsp,method = "kendall"))

[1] 0.3308049

(omega <- 2/pi*asin(tau))

[1] 0.2146404

copula2 <- tCopula(omega,dim=2,dispstr = "un",df = 2)
d1 <- cbind(fGarch::pstd(ibm, mean = est.ibm[1], sd = est.ibm[2], nu = est.ibm[3]),
            fGarch::pstd(crsp, mean = est.crsp[1], sd = est.crsp[2], nu = est.crsp[3]))
(fit1 <- fitCopula(copula2, method = "ml", optim.method = "L-BFGS-B", data = d1, start = c(omega,5), lower = c(0,2.5), upper = c(0.5,15)))

Call: fitCopula(copula2, data = d1, ... = pairlist(method = "ml", optim.method = "L-BFGS-B", start = c(omega, 
    5), lower = c(0, 2.5), upper = c(0.5, 15)))
Fit based on "maximum likelihood" and 2528 2-dimensional observations.
Copula: tCopula 
 rho.1     df 
0.4937 9.8537 
The maximized loglikelihood is 362 
Optimization converged

#Ajusta copula normal
fnorm <- fitCopula(data = d1, copula = normalCopula(-0.3, dim = 2),
method = "ml", optim.method = "BFGS", start = 0.5)
#Ajusta Gumbel
fgumbel <- fitCopula(data=d1,copula=gumbelCopula(3,dim=2),
method="ml",optim.method="BFGS",start=2)
#Ajusta Frank
ffrank <- fitCopula(data=d1,copula=frankCopula(3,dim=2),
method="ml",optim.method="BFGS",start=2)
#Ajusta Clayton
fclayton <- fitCopula(data=d1,copula=claytonCopula(3,dim=2),
method="ml",optim.method="BFGS",start=2)

Comparación gráfica

u <- d1
dem <- pempiricalCopula(u[,1],u[,2])
par(mfrow = c(3,2), mar = c(2,2,2,2))
contour(dem$x, dem$y, dem$z, main = "Cópula Empírica")
contour(tCopula(fit1@estimate[1], 
                df = round(fit1@estimate[2],0)),
        pCopula, main = "Cópula t")
contour(normalCopula(fnorm@estimate), pCopula, 
        main = "Cópula Normal")
contour(gumbelCopula(fgumbel@estimate), pCopula, 
        main = "Cópula Gumbel")
contour(frankCopula(ffrank@estimate), pCopula,
        main = "Cópula Frank")
contour(claytonCopula(fclayton@estimate), pCopula, 
        main = "Cópula Clayton")

par(mfrow=c(3,2), mar=c(2,2,2,2))
contour(kde2d(u[,1],u[,2]), main = "KDE")
contour(tCopula(fit1@estimate[1], df = fit1@estimate[2]), dCopula, 
        main = "Cópula t", 
        nlevels = 25)
contour(normalCopula(fnorm@estimate), dCopula, 
        main = "Cópula Normal", nlevels = 25)
contour(gumbelCopula(fgumbel@estimate), dCopula, 
        main = "Cópula Gumbel", nlevels = 25)
contour(frankCopula(ffrank@estimate), dCopula, 
        main = "Cópula Frank", nlevels = 25)
contour(claytonCopula(fclayton@estimate), dCopula, 
        main = "Cópula Clayton", nlevels = 25)

# AIC
# AIC Normal
2*length(fnorm@estimate) - 2*fnorm@loglik

[1] -692.3688

# AIC Gumbel
2*length(fgumbel@estimate) - 2*fgumbel@loglik

[1] -624.4514

# AIC frank
2*length(ffrank@estimate) - 2*ffrank@loglik

[1] -648.5734

# AIC Clayton
2*length(fclayton@estimate) - 2*fclayton@loglik

[1] -584.2204

# AIC t
2*length(fit1@estimate) - 2*fit1@loglik

[1] -719.9693

Dependencia de colas

Ejemplo 1

knitr::opts_chunk$set(echo = TRUE, fig.align = "center", comment = NA)
library(evd)
library(CombMSC)
library(readxl)
library(copula)
library(MASS)

datos <- read_xlsx("data/Series_Bursas_IFNB.xlsx", skip = 2) # datos originales
resCR <- read_xlsx("data/Series Bursas IFNB output.xlsx",
sheet = "Hoja2",range = "S4:BL507", col_names = T)
rescr <- resCR[,seq(2,46,by=3)]
resAC <- read_xlsx("data/Series Bursas IFNB output.xlsx",
sheet = "Hoja1",range = "V3:Bo506", col_names = T)
resac <- resAC[,seq(2,46,by=3)]
names(datos)[2] <- "Alpha_Credit"
names(datos)[3] <- "Credito_Real"

# Estimaciones empíricas
Lemp <- function(z) sum((U<=z) & (V<=z))/sum(U<=z)
Remp <- function(z) sum((U>=1-z) & (V>=1-z))/sum(U>=1-z)
L2emp <- function(z) 2*log(mean(U <= z))/log(mean((U <= z) & (V <= z))) - 1
R2emp <- function(z) 2*log(mean(U >= 1-z))/log(mean((U >= 1-z) & (V >= 1-z))) - 1

par(mfrow = c(2,2))
        n <- 5000
        x <- rnorm(n);y <- rnorm(n)
        
        # Aplicando el ejercicio:
        U <- rank(x)/(n+1)
        V <- rank(y)/(n+1)
        u <- seq(0.001,0.5,by = 0.001)
        L <- Vectorize(Lemp)(u)
        R <- Vectorize(Remp)(rev(u))
        
        
        plot(c(u, u + 0.5 - u[1]), c(L, R), type = "l", ylim = c(0,1), 
        xlab = "Cola Inferior L           Cola Superior R",
        ylab = expression(lambda))
        abline(v = 0.5, col = "grey")
        
        
        plot(U,V, main = "Rangos normalizados\n (empates = media)", pch = 16, cex = 0.5)
        plot(x, y, main = "Datos originales", pch = 16, cex = 0.5)
        plot(1:length(L), L, ylim = c(-1.2, 1.2), main = "Índices", ylab = "", xlab = "", type = "l", col = "green")
        legend("bottomright", legend = c("Superioir", "Inferior"), lty = c(1,1), col = c("red", "green"))
        lines(R, col = "red", lwd = 2)

        par(mfrow = c(2,2))
        n <- 5000
        X <- mvrnorm(n, mu = c(0,0), Sigma = matrix(c(1,0.8,0.8,1), nrow = 2))
        x <- X[,1]
        y <- X[,2]
        # Aplicando el ejercicio:
        U <- rank(x)/(n + 1)
        V <- rank(y)/(n + 1)
        u <- seq(0.001, 0.5, by = 0.001)
        L <- Vectorize(Lemp)(u)
        R <- Vectorize(Remp)(rev(u))
        
        plot(c(u, u + 0.5 - u[1]), c(L, R), type = "l", ylim = c(0,1), 
        xlab = "Cola Inferior L           Cola Superior R",
        ylab = expression(lambda))
        abline(v = 0.5, col = "grey")
        
        plot(U, V, main = "Rangos normalizados\n (empates = media)", pch = 16, cex = 0.5)
        plot(x, y, main = "Datos originales", pch = 16, cex = 0.5)
        
        plot(1:length(L), L, ylim = c(-1.2,1.2), main = "Índices", ylab = "", xlab = "", type = "l", col = "green")
        legend("bottomright", legend = c("Superioir", "Inferior"), lty = c(1,1), col = c("red", "green"))
        lines(R, col = "red", lwd = 2)

        par(mfrow = c(2,2))
        n <- 5000
        x <- rnorm(n); y <- rnorm(n)
        
        
        # Aplicando el ejercicio:
        U <- rank(x)/(n+1); V <- rank(y)/(n+1)
        u <- seq(0.001,0.5,by = 0.001)
        L <- Vectorize(L2emp)(u)
        R <- Vectorize(R2emp)(rev(u))
        
        
        plot(c(u, u + 0.5 - u[1]), c(L, R), type = "l", ylim = c(-1,1), 
        xlab = "Cola Inferior L           Cola Superior R",
        ylab = expression(xi))
        abline(v = 0.5, col = "grey")
        
        
        plot(U, V, main = "Rangos normalizados\n (empates = media)", pch = 16, cex = 0.5)
        plot(x, y, main = "Datos originales", pch = 16, cex = 0.5)
        plot(1:length(L), L, ylim = c(-1.2, 1.2), main = "Índices", ylab = "", xlab = "", type = "l", col = "green")
        legend("bottomright", legend = c("Superioir", "Inferior"), lty = c(1,1), col = c("red", "green"))
        lines(R, col = "red", lwd = 2)

par(mfrow = c(2,2))
n <- 5000
X <- mvrnorm(n, mu = c(0,0), Sigma = matrix(c(1,0.9,0.9,1), nrow = 2))
x <- X[,1]; y <- X[,2]
        
        
# Aplicando el ejercicio:
U <- rank(x)/(n + 1); V <- rank(y)/(n + 1)
u <- seq(0.001, 0.5, by = 0.001)
L <- Vectorize(L2emp)(u)
R <- Vectorize(R2emp)(rev(u))
        
plot(c(u, u + 0.5 - u[1]), c(L, R), type = "l", ylim = c(0,1), 
        xlab = "Cola Inferior L           Cola Superior R",
        ylab = expression(xi))
        abline(v = 0.5,col = "grey")
        
plot(U, V, main = "Rangos normalizados\n (empates = media)", pch = 16, cex = 0.5)
plot(x, y, main = "Datos originales", pch = 16, cex = 0.5)
plot(1:length(L), L, ylim = c(-1.2,1.2), main = "Índices", ylab = "", xlab = "", type = "l", col = "green")
legend("bottomright", legend = c("Superioir", "Inferior"), lty = c(1,1), col = c("red", "green"))
lines(R, col = "red", lwd = 2)

datos3 <- read_xlsx("data/ByA.xlsx", sheet = "Spread", range = "V8:Z638",
col_names = c("Alpha_Credit", "Crédito_Real", "Unifin", "Fin.Indep", "Mexarend"))

pares <- subsets(5, 2, 1:5)
Rl <- R2 <- list(NULL)
Ll <- L2 <- list(NULL)

par(mfrow = c(2,2))

for(i in 1:nrow(pares)){
        U <- rank(datos3[ , pares[i,1]])/(nrow(datos3)+1)
        V <- rank(datos3[ , pares[i,2]])/(nrow(datos3)+1)
        u <- seq(0.001, 0.5, by = 0.001)
        L <- Vectorize(Lemp)(u)
        R <- Vectorize(Remp)(rev(u))
        Li2 <- Vectorize(L2emp)(u)
        Ri2 <- Vectorize(R2emp)(rev(u))

        Ll[[i]] <- L
        Rl[[i]] <- R
        L2[[i]] <- Li2
        R2[[i]] <- Ri2

        print(paste(names(datos3)[pares[i,]], pares[i,]))

        # Gráfica LR  
        plot(c(u, u + 0.5 - u[1]), c(L, R), type = "l", ylim = 0:1, 
             xlab = "Cola Inferior L          Cola Superior R",
             ylab = expression(lambda),
             main = paste(names(datos3)[pares[i,1]]," y ",names(datos3)[pares[i,2]]))
        abline(v = 0.5, col = "grey")
        legend("topleft", legend = c("Fuerte","Débil"), col = c("black","navy"),lty = c(1, 0), pch = c(NA, 16), cex=0.5)

points(c(u, u + 0.5 - u[1]), c(Li2, Ri2), pch = 16, ylim = 0:1, col = "navy")
abline(v = 0.5, col = "grey")

# Datos originales
plot(datos3[, pares[i,1:2]], main = "Datos originales")
# Datos en rangos
plot(apply(datos3[, pares[i, 1:2]], 2,
           function(x)rank(x)/length(x)), 
     main = "Rangos normalizados")
nombre <- names(datos3)[pares[i,2]]

# Serie de los índices a lo largo del tiempo
plot(1:length(R), R, ylim = c(0,1), main = "Índices", 
     type = "l", col = "red", ylab = "", xlab = "")
legend("bottomright", legend = c("Sup", "Inf"), 
       lty = c(1, 1), col = c("red", "green"), cex=0.5)
points(Ri2, col = "orange")
points(Li2, col = "yellow")
lines(L, col = "green", lwd = 2)
}

[1] "Alpha_Credit 1" "Crédito_Real 2"

[1] "Alpha_Credit 1" "Unifin 3"      

[1] "Alpha_Credit 1" "Fin.Indep 4"   

[1] "Alpha_Credit 1" "Mexarend 5"    

[1] "Crédito_Real 2" "Unifin 3"      

[1] "Crédito_Real 2" "Fin.Indep 4"   

[1] "Crédito_Real 2" "Mexarend 5"    

[1] "Unifin 3"    "Fin.Indep 4"

[1] "Unifin 3"   "Mexarend 5"

[1] "Fin.Indep 4" "Mexarend 5" 

names(Ll) <- names(Rl) <- names(L2) <- names(R2) <- paste(names(datos3)[pares[,1]],names(datos3)[pares[,2]])